library(Matrix)
library(ggplot2)
library(Seurat)
library(dplyr)
library(cowplot)
library(stringr)
library(patchwork)
library(tidyverse)
library(tidyr)
library(BiocManager)
library(multtest)
library(metap)
 

#load in matrix files
library(Matrix)
Myrfhummt_dir <- "/nfs/turbo/umms-lprasov/Myrf hum scRNAseq/10x_analysis_7805-MB/Sample_7805-MB-2-GEX_CCGGCAACTG-TGTTAAACCG/filtered_feature_bc_matrix/"
Myrfhummt_barcode.path <- paste0(Myrfhummt_dir, "barcodes.tsv.gz")
Myrfhummtfeatures.path <- paste0(Myrfhummt_dir, "features.tsv.gz")
Myrfhummtmatrix.path <- paste0(Myrfhummt_dir, "matrix.mtx.gz")
Myrfhummtmat <- readMM(file = Myrfhummtmatrix.path)
Myrfhummt_feature.names = read.delim(Myrfhummtfeatures.path , header = FALSE, stringsAsFactors = FALSE)
Myrfhummt_barcode.names = read.delim(Myrfhummt_barcode.path, header = FALSE, stringsAsFactors = FALSE)
colnames(Myrfhummtmat) = Myrfhummt_barcode.names$V1
Myrfhummt_feature.names$V2 = make.unique(Myrfhummt_feature.names$V2, sep='-')
rownames(Myrfhummtmat) = Myrfhummt_feature.names$V2

#convert raw data into seurat object
Myrfhummt.seurat <- CreateSeuratObject(counts = Myrfhummtmat, project = "Myrf_mt", min.cells = 5)

Myrfhummt.seurat$Geno <- "Myrf Mutant"
Myrfhummt.seurat <- subset(Myrfhummt.seurat, subset = nFeature_RNA > 500)
Myrfhummt.seurat <- NormalizeData(Myrfhummt.seurat, verbose = FALSE)
Myrfhummt.seurat <- FindVariableFeatures(Myrfhummt.seurat, selection.method = "vst", nfeatures = 2000)

Myrfhummt.seurat[["percent.mt"]] <- PercentageFeatureSet(Myrfhummt.seurat, pattern = "^mt-")

VlnPlot(Myrfhummt.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
Myrfhummtselect.seurat <- subset(Myrfhummt.seurat, subset = nFeature_RNA > 200 & percent.mt < 15)
VlnPlot(Myrfhummtselect.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot1_mut <- FeatureScatter(Myrfhummtselect.seurat, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2_mut <- FeatureScatter(Myrfhummtselect.seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1_mut, plot2_mut))

#combine the seurat objects
#Dim default is 30, maybe 20 is more stringent... ask Leonard
Myrfhumselect.anchors <- FindIntegrationAnchors(object.list = list(Myrfhumwt.seurat, Myrfhummtselect.seurat), dims = 1:20)

Myrfhumgroup.seurat <- IntegrateData(anchorset = Myrfhumselect.anchors, dims = 1:20)

#The object Myrfgroup is the one that I use most frequently.  It has the ability to sort by genotype or by group (WT vs MUT).  It also has the clusters organized such that all the RPC are together, etc.
DefaultAssay(Myrfhumgroup.seurat) <- "integrated"

# Run the standard workflow for visualization and clustering
Myrfhumgroup.seurat  <- ScaleData(Myrfhumgroup.seurat, verbose = FALSE)
Myrfhumgroup.seurat  <- RunPCA(Myrfhumgroup.seurat, npcs = 30, verbose = FALSE)
# UMAP and Clustering
Myrfhumgroup.seurat  <- RunUMAP(Myrfhumgroup.seurat, reduction = "pca", dims = 1:20)
Myrfhumgroup.seurat  <- FindNeighbors(Myrfhumgroup.seurat, reduction = "pca", dims = 1:20)
Myrfhumgroup.seurat <- FindClusters(Myrfhumgroup.seurat, resolution = 0.5)

#view clusters
DimPlot(Myrfhumgroup.seurat, reduction = "umap", label = TRUE)

#check gene expression of MITF and MYRF
FeaturePlot(Myrfhumgroup.seurat, features = c("Mitf"))
FeaturePlot(Myrfhumgroup.seurat, features = c("Mitf", "Myrf"))

#determining clusters
Myrfhummarkers <- FindAllMarkers(Myrfhumgroup.seurat, only.pos = TRUE)

##Anchoring datasets
library(Seurat)
library(dplyr)

#cowplot is an add on package to ggplot to make publication quality images
library(cowplot)
#load in matrix files
library(Matrix)
e13WT_matrix_dir = "/nfs/turbo/umms-lprasov/Prasov/Run_2961/Sample_Myrf_wt/outs/filtered_feature_bc_matrix/"
e13WT_barcode.path <- paste0(e13WT_matrix_dir, "barcodes.tsv.gz")
e13WT_features.path <- paste0(e13WT_matrix_dir, "features.tsv.gz")
e13WT_matrix.path <- paste0(e13WT_matrix_dir, "matrix.mtx.gz")
e13WT_mat <- readMM(file = e13WT_matrix.path)
e13WT_feature.names = read.delim(e13WT_features.path, header = FALSE, stringsAsFactors = FALSE)
e13WT_barcode.names = read.delim(e13WT_barcode.path, header = FALSE, stringsAsFactors = FALSE)
colnames(e13WT_mat) = e13WT_barcode.names$V1
e13WT_feature.names$V2 = make.unique(e13WT_feature.names$V2, sep='-')
rownames(e13WT_mat) = e13WT_feature.names$V2
load("/nfs/turbo/umms-lprasov/MyrfhumSeurat/.RData")
e15WT_matrix_dir = "/nfs/turbo/umms-lprasov/NovaA161_e15scseq/Sample_Myrf_fl_fl/outs/filtered_feature_bc_matrix/"
e15WT_barcode.path <- paste0(e15WT_matrix_dir, "barcodes.tsv.gz")
e15WT_features.path <- paste0(e15WT_matrix_dir, "features.tsv.gz")
e15WT_matrix.path <- paste0(e15WT_matrix_dir, "matrix.mtx.gz")
e15WT_mat <- readMM(file = e15WT_matrix.path)
e15WT_feature.names = read.delim(e15WT_features.path, header = FALSE, stringsAsFactors = FALSE)
e15WT_barcode.names = read.delim(e15WT_barcode.path, header = FALSE, stringsAsFactors = FALSE)
colnames(e15WT_mat) = e15WT_barcode.names$V1
e15WT_feature.names$V2 = make.unique(e15WT_feature.names$V2, sep='-')
rownames(e15WT_mat) = e15WT_feature.names$V2

#convert raw data into seurat object
e15WT.seurat <- CreateSeuratObject(counts = e15WT_mat, project = "e15_Myrf_wt", min.cells = 5)
e15MUT_matrix_dir = "/nfs/turbo/umms-lprasov/NovaA161_e15scseq/Sample_Rxcre_Myrf_fl_fl/outs/filtered_feature_bc_matrix/"
e15MUT_barcode.path <- paste0(e15MUT_matrix_dir, "barcodes.tsv.gz")
e15MUT_features.path <- paste0(e15MUT_matrix_dir, "features.tsv.gz")
e15MUT_matrix.path <- paste0(e15MUT_matrix_dir, "matrix.mtx.gz")
e15MUT_mat <- readMM(file = e15MUT_matrix.path)
e15MUT_feature.names = read.delim(e15MUT_features.path, header = FALSE, stringsAsFactors = FALSE)
e15MUT_barcode.names = read.delim(e15MUT_barcode.path, header = FALSE, stringsAsFactors = FALSE)
colnames(e15MUT_mat) = e15MUT_barcode.names$V1
e15MUT_feature.names$V2 = make.unique(e15MUT_feature.names$V2, sep='-')
rownames(e15MUT_mat) = e15MUT_feature.names$V2

#convert raw data into seurat object
e15MUT.seurat <- CreateSeuratObject(counts = e15MUT_mat, project = "e15_Rxcre_Myrf_mut", min.cells = 5)
P0WT_matrix_dir = "/nfs/turbo/umms-lprasov/Nova_170/NovaA_170/Sample_Myrf_wt/outs/filtered_feature_bc_matrix/"
P0WT_barcode.path <- paste0(P0WT_matrix_dir, "barcodes.tsv.gz")
P0WT_features.path <- paste0(P0WT_matrix_dir, "features.tsv.gz")
P0WT_matrix.path <- paste0(P0WT_matrix_dir, "matrix.mtx.gz")
P0WT_mat <- readMM(file = P0WT_matrix.path)
P0WT_feature.names = read.delim(P0WT_features.path, header = FALSE, stringsAsFactors = FALSE)
P0WT_barcode.names = read.delim(P0WT_barcode.path, header = FALSE, stringsAsFactors = FALSE)
colnames(P0WT_mat) = P0WT_barcode.names$V1
P0WT_feature.names$V2 = make.unique(P0WT_feature.names$V2, sep='-')
rownames(P0WT_mat) = P0WT_feature.names$V2

#convert raw data into seurat object
P0WT.seurat <- CreateSeuratObject(counts = P0WT_mat, project = "P0_Myrf_wt", min.cells = 5)
P0MUT_matrix_dir = "/nfs/turbo/umms-lprasov/Nova_170/NovaA_170/Sample_Rxcre_Myrf_mut/outs/filtered_feature_bc_matrix/"
P0MUT_barcode.path <- paste0(P0MUT_matrix_dir, "barcodes.tsv.gz")
P0MUT_features.path <- paste0(P0MUT_matrix_dir, "features.tsv.gz")
P0MUT_matrix.path <- paste0(P0MUT_matrix_dir, "matrix.mtx.gz")
P0MUT_mat <- readMM(file = P0MUT_matrix.path)
P0MUT_feature.names = read.delim(P0MUT_features.path, header = FALSE, stringsAsFactors = FALSE)
P0MUT_barcode.names = read.delim(P0MUT_barcode.path, header = FALSE, stringsAsFactors = FALSE)
colnames(P0MUT_mat) = P0MUT_barcode.names$V1
P0MUT_feature.names$V2 = make.unique(P0MUT_feature.names$V2, sep='-')
rownames(P0MUT_mat) = P0MUT_feature.names$V2

#convert raw data into seurat object
P0MUT.seurat <- CreateSeuratObject(counts = P0MUT_mat, project = "P0_Rxcre_Myrf_mut", min.cells = 5)
e15WT.seurat$Geno <- "e15_WT"
e15WT.seurat <- subset(e15WT.seurat, subset = nFeature_RNA > 500)
e15WT.seurat <- NormalizeData(e15WT.seurat, verbose = FALSE)
e15WT.seurat <- FindVariableFeatures(e15WT.seurat, selection.method = "vst", nfeatures = 2000)
e15MUT.seurat$Geno <- "e15_MUT"
e15MUT.seurat <- subset(e15MUT.seurat, subset = nFeature_RNA > 500)
e15MUT.seurat <- NormalizeData(e15MUT.seurat, verbose = FALSE)
e15MUT.seurat <- FindVariableFeatures(e15MUT.seurat, selection.method = "vst", nfeatures = 2000)

P0WT.seurat$Geno <- "P0_WT"
P0WT.seurat <- subset(P0WT.seurat, subset = nFeature_RNA > 500)
P0WT.seurat <- NormalizeData(P0WT.seurat, verbose = FALSE)
P0WT.seurat <- FindVariableFeatures(P0WT.seurat, selection.method = "vst", nfeatures = 2000)

P0MUT.seurat$Geno <- "P0-MUT"
P0MUT.seurat <- subset(P0MUT.seurat, subset = nFeature_RNA > 500)
P0MUT.seurat <- NormalizeData(P0MUT.seurat, verbose = FALSE)
P0MUT.seurat <- FindVariableFeatures(P0MUT.seurat, selection.method = "vst", nfeatures = 2000)

e15WT.seurat[["percent.mt"]] <- PercentageFeatureSet(e15WT.seurat, pattern = "^mt-")
e15MUT.seurat[["percent.mt"]] <- PercentageFeatureSet(e15MUT.seurat, pattern = "^mt-")
P0WT.seurat[["percent.mt"]] <- PercentageFeatureSet(P0WT.seurat, pattern = "^mt-")
P0MUT.seurat[["percent.mt"]] <- PercentageFeatureSet(P0MUT.seurat, pattern = "^mt-")
e15WTselect.seurat <- subset(e15WT.seurat, subset = nFeature_RNA > 200 & percent.mt < 15)
e15MUTselect.seurat <- subset(e15MUT.seurat, subset = nFeature_RNA > 200 & percent.mt < 15)
P0WTselect.seurat <- subset(P0WT.seurat, subset = nFeature_RNA > 200 & percent.mt < 15)
P0MUTselect.seurat <- subset(P0MUT.seurat, subset = nFeature_RNA > 200 & percent.mt < 15)
Myrfhumwt.seurat[["percent.mt"]] <- PercentageFeatureSet(Myrfhumwt.seurat, pattern = "^mt-")
Myrfhumwtselect.seurat <- subset(Myrfhumwt.seurat, subset = nFeature_RNA > 200 & percent.mt < 15)

#Anchor the datasets together
Myrfhum.anchors <- FindIntegrationAnchors(object.list = list(e15WTselect.seurat, e15MUTselect.seurat, P0WTselect.seurat, P0MUTselect.seurat, Myrfhumwtselect.seurat, Myrfhummtselect.seurat), dims = 1:20)
?save.image()

#Integrate datasets
Myrfhumgroup.seurat <- IntegrateData(anchorset = Myrfhum.anchors, dims = 1:20)


# Run the standard workflow for visualization and clustering
Myrfhumgroup.seurat  <- ScaleData(Myrfhumgroup.seurat, verbose = FALSE)
Myrfhumgroup.seurat  <- RunPCA(Myrfhumgroup.seurat, npcs = 30, verbose = FALSE)

# UMAP and Clustering
Myrfhumgroup.seurat  <- RunUMAP(Myrfhumgroup.seurat, reduction = "pca", dims = 1:20)
Myrfhumgroup.seurat  <- FindNeighbors(Myrfhumgroup.seurat, reduction = "pca", dims = 1:20)
Myrfhumgroup.seurat <- FindClusters(Myrfhumgroup.seurat, resolution = 0.5)

#Visualize clusters
p1 <- DimPlot(Myrfhumgroup.seurat, reduction = "umap", group.by = "Geno")
p2 <- DimPlot(Myrfhumgroup.seurat, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

DimPlot(Myrfhumgroup.seurat, reduction = "umap", split.by = "Geno")

## Visualizing and looking into potential causes of the 0 cluster
FeaturePlot(Myrfhumgroup.seurat, 'nCount_RNA', order=T)
FeaturePlot(Myrfhumgroup.seurat, 'nFeature_RNA', order=T)
FeaturePlot(Myrfhumgroup.seurat, 'percent.mt', order=T)
VlnPlot(Myrfhumgroup.seurat, 'nFeature_RNA')
VlnPlot(Myrfhumgroup.seurat, 'nFeature_RNA', pt.size = 0)&NoLegend()
DimPlot(Myrfhumgroup.seurat, label = T)&NoLegend()
VlnPlot(Myrfhumgroup.seurat, 'nCount_RNA', pt.size = 0)&NoLegend()
VlnPlot(Myrfhumgroup.seurat, 'nCount_RNA', pt.size = 0, ymax = 8000)&NoLegend()

##Ran this FindAllMarkers on the non-subsetted group
# For performing differential expression after integration, we switch back to the original data
## FindAllMarkers --> better for just one condition or one group
DefaultAssay(Myrfhumgroup.seurat) <- "RNA"
## Default statistical test here is wilcox but I have seen other people use DESeq2 instead...
all.selected.markers <- FindAllMarkers (Myrfhumgroup.seurat, only.pos = TRUE)

##FindConservedMarkers
DefaultAssay(Myrfhumgroup.seurat) <- "RNA"
markers_cluster0 <- FindConservedMarkers(Myrfhumgroup.seurat, ident.1 = 0, grouping.var = "Geno", verbose = FALSE)

#FindConservedMarkers on subset data
DefaultAssay(Myrfhumgroup_sub.seurat) <- "RNA"
markers_cluster0 <- FindConservedMarkers(Myrfhumgroup_sub.seurat, ident.1 = "0", grouping.var = "Geno", verbose = FALSE)
markers_cluster1 <- FindConservedMarkers(Myrfhumgroup_sub.seurat, ident.1 = "1", grouping.var = "Geno", verbose = FALSE)

AllClus <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25)

for (i in AllClus) {
  DefaultAssay(Myrfhumgroup_sub.seurat) <- "RNA"
  gmarker.markers <- FindConservedMarkers(Myrfhumgroup_sub.seurat, ident.1 = i, grouping.var = "Geno", verbose = FALSE)
  assign(paste0("g", i, ".markers"), gmarker.markers)
}
#cleanup environment
rm("gmarker.markers", "AllClus")

save(gmarker.markers, file = "gmarker_markers.csv")

## To save FindAllMarkers as a excel file
write.csv(all.selected.markers, "Myrfhum_allselectedmarkers.csv")
head(all.selected.markers)

##Trying to match up marker genes to clusters
p1 <- DimPlot(Myrfhumgroup.seurat, reduction = "umap", group.by = "Geno")
p2 <- DimPlot(Myrfhumgroup.seurat, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

FeaturePlot(Myrfhumgroup.seurat, features = c("Rho", "Gnat1", "Cnga1", "Nrl", "Nr2e3"), min.cutoff = "q9")

FeaturePlot(Myrfhumgroup.seurat, features = c("Sncg", "Pou4f1", "Nrn1", "Slc17a6"), min.cutoff = "q9")

FeaturePlot(Myrfhumgroup.seurat, features = c("Lhx1", "Calb1", "Gja10", "Onecut1"), min.cutoff = "q9")

## Name clusters based on marker genes
# Myrfhumgroup.seurat<- RenameIdents(Myrfgroup.seurat, '0' = "??", '1' = "MG", '2' ="", '3' = "Rods", '4' = "Ama4", '5' = "PhRe5", '6' = "Neu6", '7' = "POM7", '8' = "RGC", '9' = "PhRec9", '10' = "RPC10", '11' = "HC"", '12' = "ONH12", '13' = "RGC13", '14' = "RPEall", '15' = "RPC15", '16' = "POM16", '17' = "POM17", '18' = "RPC18", '19' = "BV19", '20' = "???", '21' = "???", '22' = "Cones", '23' = "???", '24' = "???", '25' = "???")

## Remove cluster 0
# To subset and remove single cluster and keep the remaining clusters for new analysis
Myrfhumgroup_sub.seurat <- subset(Myrfhumgroup.seurat, idents = 0, invert = TRUE)

##Redoing clustering and analysis on subset group without cluster 0
#The object Myrfgroup is the one that I use most frequently.  It has the ability to sort by genotype or by group (WT vs MUT).  It also has the clusters organized such that all the RPC are together, etc.
DefaultAssay(Myrfhumgroup_sub.seurat) <- "integrated"

# Run the standard workflow for visualization and clustering
Myrfhumgroup_sub.seurat  <- ScaleData(Myrfhumgroup_sub.seurat, verbose = FALSE)
Myrfhumgroup_sub.seurat  <- RunPCA(Myrfhumgroup_sub.seurat, npcs = 30, verbose = FALSE)
# UMAP and Clustering
Myrfhumgroup_sub.seurat  <- RunUMAP(Myrfhumgroup_sub.seurat, reduction = "pca", dims = 1:20)
Myrfhumgroup_sub.seurat  <- FindNeighbors(Myrfhumgroup_sub.seurat, reduction = "pca", dims = 1:20)
Myrfhumgroup_sub.seurat <- FindClusters(Myrfhumgroup_sub.seurat, resolution = 0.5)

#visualize clsuters
p1 <- DimPlot(Myrfhumgroup_sub.seurat, reduction = "umap", group.by = "Geno")
p2 <- DimPlot(Myrfhumgroup_sub.seurat, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

DimPlot(Myrfhumgroup_sub.seurat, reduction = "umap", split.by = "Geno")

#Save subset dataset
save(Myrfhumgroup_sub.seurat, file = "Myrfhumgroup_sub.seurat.Rdata")
save(Cluster.merged.seurat, file = "Cluster.merged.seurat.Rdata")

#To get number of cells in clusters not split by geno, this will list give the cluster order
numcells_clust <- table(Cluster.merged.seurat@active.ident)
write.csv(numcells_clust, "numcells_clust.csv")

#To get total number cells per geno
numcells_bygeno <- table(Cluster.merged.seurat$Geno)
write.csv(numcells_bygeno, "numcells_bygeno.csv")

#To get total number of cells per cluster split by geno
numcells_perclust_bygeno <- table(Cluster.merged.seurat@active.ident, Cluster.merged.seurat$Geno)
write.csv(numcells_perclust_bygeno, "TotalCells_perClust_bygeno.csv")

### Ran on 7/10/23 evening @ 10:45pm
## FindAllMarkers on subset dataset --> better for just one condition or one group
DefaultAssay(Myrfhumgroup_sub.seurat) <- "RNA"
## Default statistical test here is wilcox but I have seen other people use DESeq2 instead...
all.subselected.markers <- FindAllMarkers(Myrfhumgroup_sub.seurat, only.pos = TRUE)
## To save FindAllMarkers as a excel file
write.csv(all.subselected.markers, "Myrfhum_sub_allselectedmarkers.csv")


#Checking common quality controls
FeaturePlot(Myrfhumgroup_sub.seurat, 'nCount_RNA', order=T)
FeaturePlot(Myrfhumgroup_sub.seurat, 'nFeature_RNA', order=T)
FeaturePlot(Myrfhumgroup_sub.seurat, 'percent.mt', order=T)
VlnPlot(Myrfhumgroup_sub.seurat, 'nFeature_RNA')
VlnPlot(Myrfhumgroup_sub.seurat, 'nFeature_RNA', pt.size = 0)&NoLegend()
DimPlot(Myrfhumgroup_sub.seurat, label = T)&NoLegend()
VlnPlot(Myrfhumgroup_sub.seurat, 'nCount_RNA', pt.size = 0)&NoLegend()
VlnPlot(Myrfhumgroup_sub.seurat, 'nCount_RNA', pt.size = 0, ymax = 8000)&NoLegend()

#Labeling cell clusters
#ama
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Ptf1a", "Onecut2", "Prdm13"), min.cutoff = "q9")
#RBC
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Hbb-y", "Hba-a2"), min.cutoff = "q9")
#BV
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Egfl7", "Ramp2", "Igfbp7"), min.cutoff = "q9")
#RPE
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Dct", "Mitf", "Pmel", "Tyr", "Ttr"), min.cutoff = "q9")
#Mel
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Dct", "Mitf", "Pmel", "Tyr", "Pax3", "Sox10", "Foxd3", "Ptgds"), min.cutoff = "q9")
#RPC
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Vsx2", "Lhx2", "Sorbs2", "Eno1", "Ccnd1", "Fgf15"), min.cutoff = "q9")
#POM
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Mgp", "Pitx2", "Irx1"), min.cutoff = "q9")
#LSC
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Klf5", "Krt5", "Krt7", "S100a6"), min.cutoff = "q9")
#WBC
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("C1qc", "C1qb"), min.cutoff = "q9")
#Cor
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Osr1", "Klf4", "Foxl2", "Dlk1", "Lgals1", "Mfap5", "Eln", "Postn"), min.cutoff = "q9")
#EOM
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Pax7", "Myog", "Mstn", "Actc1", "Cdkn1c"), min.cutoff = "q9")
#CilEp
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Kera", "Matn4", "Aqp1", "Ald3a1"), min.cutoff = "q9")
#Peri
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Ecm1", "Rgs5", "Tbx18"), min.cutoff = "q9")
#PhRep
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Lhx4", "Otx2", "Crx", "Prdm1"), min.cutoff = "q9")
#Neur
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Foxn4", "Notch1", "Rbpj"), min.cutoff = "q9")
#Conj
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Postn", "Col3a1", "Tnmd", "Col1a1", "Col1a2"), min.cutoff = "q9")
#RGC
FeaturePlot(Myrfhumgroup_sub.seurat, features = c("Isl1", "Pou4f2", "Atho7", "Tubb3", "Stmn2", "Gap43", "Sox11", "Nnat"), min.cutoff = "q9")
Myrfhumgroup_sub.seurat<- RenameIdents(Myrfhumgroup_sub.seurat, '5' = "RGC5")
Myrfhumgroup_sub.seurat<- RenameIdents(Myrfhumgroup_sub.seurat, 'PhRep21' = "RGC21")
Myrfhumgroup_sub.seurat<- RenameIdents(Myrfhumgroup_sub.seurat, 'RGC6' = "COR6")

#Visualize with new cluster names
p1 <- DimPlot(Myrfhumgroup_sub.seurat, reduction = "umap", group.by = "Geno")
p2 <- DimPlot(Myrfhumgroup_sub.seurat, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
#Change cluster names
Myrfhumgroup_sub.seurat<- RenameIdents(Myrfhumgroup_sub.seurat, '0' = "RPC0", '1' = "PhRep1", '2' = "RPC2", '3' = "RPC3", '4' = "NEUR4", '5' = "RGC5", '6' = 'COR6', '7' = "COR7", '8' = "NEUR8", '9' = "POM9", '10' = "RPE10", '11' = "AMA11", '12' = "RPC12", '13' = "POM13", '14' = "BV14", '15' = "CilEp15", '16' = "WBC16", '17' = "17", '18' = "EOM18", '19' = "RBC19", '20' = "LSC20", '21' = "RGC21", '22' = "MEL22", '23' = "23", '24' = "24", '25' = "CilEp25")
Myrfhumgroup_sub.seurat <- RenameIdents(Myrfhumgroup_sub.seurat, 'POM9' = "PERI9", '17' = "CONJ17")

#To merge clusters
Myrfhumgroup_sub.seurat <- RenameIdents(Myrfhumgroup_sub.seurat, 'RPC2' = "RPC0", 'RPC3' = "RPC0", 'RPC12'= "RPC0", 'COR7' = "COR6", 'CilEp25' = "CilEp15", 'NEUR8' = "NEUR4")
#To confirm merge, look at umap
DimPlot(object = Myrfhumgroup_sub.seurat, pt.size = 0.5, label = T, reduction = 'umap')
DimPlot(RPC.merged, reduction = "umap", split.by = "Geno")
#Then I renamed as RPC instead of RPC0
Myrfhumgroup_sub.seurat <- RenameIdents(object = Myrfhumgroup_sub.seurat, 'RPC0' = "RPC", 'COR6' = "COR", 'CilEp15' = "CilEp", 'NEUR4' = "NEUR", 'PhRep1' = "PhRep", 'RGC5' = "RGC", 'PERI9' = "PERI", 'RPE10' = "RPE", 'AMA11' = "AMA", 'POM13' = "POM", 'BV14' = "BV", 'WBC16'= "WBC", 'CONJ17' = "CONJ", 'EOM18' = "EOM", 'RBC19' = "RBC", 'LSC20' = "LSC", 'RGC21' = "AMA", 'MEL22' = "MEL", '23' = "Immature Oligo", '24' = "ONH")

DimPlot(Myrfhumgroup_sub.seurat, reduction = "umap", split.by = "Geno")

Cluster.merged.seurat <- RenameIdents(object = Cluster.merged.seurat, '21' = "AMA", '24' = "ONH", '23' = "Oligos")

#Cell Cycle analysis
#Cell cycle genes in Seurat are stored as all upper case, the 10x output is Sentence case(cap first letter), you have to convert the list
library(stringr)
m.s.genes <- str_to_title(cc.genes.updated.2019$s.genes)
m.g2m.genes <- str_to_title(cc.genes.updated.2019$g2m.genes)
hMYRFcellcycle <- CellCycleScoring(Myrfhumgroup_sub.seurat, s.features = m.s.genes, g2m.features = m.g2m.genes, set.ident = TRUE)
hMYRFcellcycle <- ScaleData(hMYRFcellcycle, verbose = FALSE)
hMYRFcellcycle <- RunPCA(hMYRFcellcycle, features = c(m.s.genes, m.g2m.genes))
DimPlot(hMYRFcellcycle, reduction = "umap", split.by = "Geno")


#Trying to visualize plots for genes by cell type
VlnPlot(Myrfhumgroup_sub.seurat, "Myrf", pt.size = 0, group.by = "Geno")
plots <- VlnPlot(Myrfhumgroup_sub.seurat, features = c("Myrf"), split.by = "Geno", group.by = "Idents", pt.size = 0, combine = FALSE)
wrap_plots(plots = plots, ncol = 1)

#Creating subsets for clusters with unknown identities
# To subset and remove single cluster and keep the remaining clusters for new analysis
# RPE subset
Myrfhumgroup_RPE.seurat <- subset(Myrfhumgroup_sub.seurat, idents = 10, invert = FALSE)
Myrfhumgroup_RPE.seurat <- RenameIdents(Myrfhumgroup_RPE.seurat, '10' = "RPE")
DefaultAssay(Myrfhumgroup_RPE.seurat) <- "RNA"
Myrfhumgroup_RPE.seurat$celltype.geno <- paste(Idents(Myrfhumgroup_RPE.seurat), Myrfhumgroup_RPE.seurat$Geno, sep = "_")
Myrfhumgroup_RPE.seurat$celltype <- Idents(Myrfhumgroup_RPE.seurat)
Idents(Myrfhumgroup_RPE.seurat) <- "celltype.geno"
#Looking at top markers in WT vs. MUT Myrf C-term mice
myrfloss.RPE <- FindMarkers(Myrfhumgroup_RPE.seurat, ident.1 = "RPE_Myrf Mutant", ident.2 = "RPE_Myrf Control", verbose = FALSE)
head(myrfloss.RPE, n = 15)


library(ggplot2)
library(cowplot)
theme_set(theme_cowplot())
Myrfhumgroup_RPE.seurat <- subset(Myrfhumgroup_RPE.seurat, idents = "RPE", invert = FALSE)
RPE.cells <- subset(Myrfhumgroup_RPE.seurat, idents = "RPE")
Idents(RPE.cells) <- "Geno"
avg.RPE.cells <- as.data.frame(log1p(AverageExpression(RPE.cells, verbose = FALSE)$RNA))
avg.RPE.cells$gene <- rownames(avg.RPE.cells)

colnames(avg.RPE.cells)[5] <- "hMyrf_WT"
colnames(avg.RPE.cells)[6] <- "hMyrf_MUT"

#upregulated genes
genes.to.label = c("Hist1h4d", "Hist1h4h", "Hist1h1e" , "Hist1h3e", "Hist1h3c", "Cdk8", "Nnat", "Tuba1a", "Hist3h2ba", "Tubb2b", "Ckb", "Hist1h1d", "Hist1h1b", "Hes5")
p1 <- ggplot(avg.RPE.cells, aes(hMyrf_WT, hMyrf_MUT)) + geom_point() + ggtitle("RPE Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p1

FeaturePlot(Myrfhumgroup_RPE.seurat, features = c("Hist1h4d", "Hist1h4h", "Hist1h1e" , "Hist1h3e", "Hist1h3c", "Cdk8", "Nnat", "Tuba1a", "Hist3h2ba", "Tubb2b", "Ckb", "Hist1h1d", "Hist1h1b", "Hes5"), split.by = "Geno", max.cutoff = 3,
            cols = c("grey", "red"))

#downregulated genes
genes.to.label = c("Ankfn1", "Pde1c", "Otx2os1", "Trpm3", "Col3a1", "Serpine3", "Col2a1", "Ttn", "Rora", "Rbms1", "Col1a1")
p2 <- ggplot(avg.RPE.cells, aes(hMyrf_WT, hMyrf_MUT)) + geom_point() + ggtitle("RPE Cells")
p2 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p2

#Looking at specific genes in RPE cluster
DimPlot(Myrfhumgroup_RPE.seurat, reduction = "umap", split.by = "Geno")
plots <- VlnPlot(Myrfhumgroup_RPE.seurat, features = c("Myrf", "Mitf", "Ermn", "Tmem98"), group.by = "Geno",
                 pt.size = 0, combine = FALSE)
plots
wrap_plots(plots = plots, ncol = 1)

FeaturePlot(Myrfhumgroup_RPE.seurat, features = "Tmem98", split.by = "Geno", max.cutoff = 3, cols = c("grey", "red"))
VlnPlot(Myrfhumgroup_RPE.seurat, features = c("Tmem98"), split.by = "celltype", group.by = "Geno",
        pt.size = 0, combine = FALSE)

up_plots <- VlnPlot(Myrfhumgroup_RPE.seurat, features = c("Hist1h4d", "Hist1h4h", "Hist1h1e" , "Hist1h3e", "Hist1h3c", "Cdk8", "Nnat", "Tuba1a", "Hist3h2ba", "Tubb2b", "Ckb", "Hist1h1d"), split.by = "celltype", group.by = "Geno",
                 pt.size = 0, combine = TRUE)
wrap_plots(plots = up_plots, ncol = 1)

down_plots <- VlnPlot(Myrfhumgroup_RPE.seurat, features = c("Ankfn1", "Pde1c", "Otx2os1", "Trpm3", "Col3a1", "Serpine3", "Col2a1", "Ttn", "Rora", "Rbms1", "Col1a1"), split.by = "celltype", group.by = "Geno",
                 pt.size = 0, combine = TRUE)
wrap_plots(plots = down_plots, ncol = 1)

#Cell cycle analysis RPE
library(stringr)
m.s.genes <- str_to_title(cc.genes.updated.2019$s.genes)
m.g2m.genes <- str_to_title(cc.genes.updated.2019$g2m.genes)
hMyrf_RPEcellcycle <- CellCycleScoring(Myrfhumgroup_RPE.seurat, s.features = m.s.genes, g2m.features = m.g2m.genes, set.ident = TRUE)
hMyrf_RPEcellcycle <- ScaleData(hMyrf_RPEcellcycle, verbose = FALSE)
hMyrf_RPEcellcycle <- RunPCA(hMyrf_RPEcellcycle, features = c(m.s.genes, m.g2m.genes))
DimPlot(hMyrf_RPEcellcycle, reduction = "umap", split.by = "Geno")
table(hMyrf_RPEcellcycle@active.ident, hMyrf_RPEcellcycle$Geno)


## Cluster RPC
Myrfhumgroup_RPC.seurat <-subset(Cluster.merged.seurat, idents = 'RPC', invert = FALSE)
DefaultAssay(Myrfhumgroup_RPC.seurat) <- "RNA"
DimPlot(Myrfhumgroup_RPC.seurat, reduction = "umap", split.by = "Geno")
Myrfhumgroup_RPC.seurat$celltype.geno <- paste(Idents(Myrfhumgroup_RPC.seurat), Myrfhumgroup_RPC.seurat$Geno, sep = "_")
Myrfhumgroup_RPC.seurat$celltype <- Idents(Myrfhumgroup_RPC.seurat)
Idents(Myrfhumgroup_RPC.seurat) <- "celltype.geno"
#Looking at top markers in WT vs. MUT Myrf C-term mice
myrfloss.RPC <- FindMarkers(Myrfhumgroup_RPC.seurat, ident.1 = "RPC_Myrf Mutant", ident.2 = "RPC_Myrf Control", verbose = FALSE)
head(myrfloss.RPC, n = 15)

write.csv(myrfloss.RPC, file = "myrfloss_RPC_topmarkers.csv")
save(Myrfhumgroup_RPC.seurat, file = "Myrfhumgroup_RPC.seurat.Rdata")

library(ggplot2)
library(cowplot)
theme_set(theme_cowplot())
RPC.cells <- subset(Myrfhumgroup_RPC.seurat, idents = "RPC")
Idents(RPC.cells) <- "Geno"
avg.RPC.cells <- as.data.frame(log1p(AverageExpression(RPC.cells, verbose = FALSE)$RNA))
avg.RPC.cells$gene <- rownames(avg.RPC.cells)

colnames(avg.RPC.cells)[5] <- "hMyrf_WT"
colnames(avg.RPC.cells)[6] <- "hMyrf_MUT"

write.csv(avg.RPC.cells, file = "avgRPCcells_DEGs.csv")

#upregulated genes
genes.to.label = c("Xist",	"Gm47283", "Hist1h3e",	"Hist1h2bn",	"Hist1h2bl",	"Cdk8",	"Hist2h2bb",	"Hist1h4d",	"Hist1h3c",	"Hist1h4h",	"Hist1h2bb",	"Hist1h2be",	"Fgf15",	"Hist1h2bg",	"Hist1h3i",	"Hist1h3f",	"Hist3h2ba",	"Hist1h1e",	"Hist1h2bm",	"Hist1h1d")
p1 <- ggplot(avg.RPE.cells, aes(hMyrf_WT, hMyrf_MUT)) + geom_point() + ggtitle("RPE Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p1

#downregulated genes
genes.to.label = c("Col3a1",	"Col1a1",	"Col1a2",	"Ttn",	"H19",	"Igf2",	"Postn",	"Camk1d",	"Uba52",	"Col2a1",	"Col11a1",	"Dcn",	"Gm10076",	"Eln",	"Ahnak",	"Acta1",	"Gm42418",	"Igfbp5",	"Gm19951",	"Mgp")
p2 <- ggplot(avg.RPE.cells, aes(hMyrf_WT, hMyrf_MUT)) + geom_point() + ggtitle("RPE Cells")
p2 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p2

#Cell Cycle Analysis in RPC
library(stringr)
m.s.genes <- str_to_title(cc.genes.updated.2019$s.genes)
m.g2m.genes <- str_to_title(cc.genes.updated.2019$g2m.genes)
hMyrf_RPCcellcycle <- CellCycleScoring(Myrfhumgroup_RPC.seurat, s.features = m.s.genes, g2m.features = m.g2m.genes, set.ident = TRUE)
hMyrf_RPCcellcycle <- ScaleData(hMyrf_RPCcellcycle, verbose = FALSE)
hMyrf_RPCcellcycle <- RunPCA(hMyrf_RPCcellcycle, features = c(m.s.genes, m.g2m.genes))
DimPlot(hMyrf_RPCcellcycle, reduction = "umap", split.by = "Geno")
table(hMyrf_RPCcellcycle@active.ident,hMyrf_RPCcellcycle$Geno)8


#Reading in all selected markers
Myrfhum_sub_allselectedmarkers.df <- read.csv("Myrfhum_sub_allselectedmarkers.csv")

#Trying to make dotplots for markers in unknown clusters
Idents(Myrfhumgroup_sub.seurat) <- factor(Idents(Myrfhumgroup_sub.seurat), levels = c("RPC0","PhRep1","RPC2","RPC3", "NEUR4", "RGC5", 'COR6', "COR7", "NEUR8", "POM9", "RPE10", "AMA11", "RPC12", "POM13", "BV14", "CilEp15", "WBC16", "17", "EOM18", "RBC19", "LSC20", "RGC21", "MEL22", "23", "24", "CilEp25"))
markers.to.plot <- c("Hba-a1", "Hbb-bt", "Hbb-y", "Ccnd1", "Dapl1", "Fgf15", "S100a6", "Krt15", "Krt5", "C1qb", "C1qa", "Pf4", "Ramp2", "Egfl7", "Plvap", "Igfbp7", "Dct", "Ptgds", "Pmel", "Ttr", "Hes6", "Neurog2", "Atoh7", "Gngt2", "Fabp7", "Neurod1", "Tubb3", "Stmn2", "Gap43", "Ptf1a", "Sox11", "Nnat", "Igfbp2", "Twist1", "Col1a1", "Col1a2", "Mgp", "Spp1", "Dlk1", "Lgals1", "Mfap5", "Eln", "Postn", "Actc1", "Acta1", "Cdkn1c")
DotPlot(Myrfhumgroup_sub.seurat, features = markers.to.plot, cols = c("blue", "red", "green", "black", "grey", "pink"), dot.scale = 8) +
  RotatedAxis()

Idents(Myrfhumgroup_sub.seurat) <- factor(Idents(Myrfhumgroup_sub.seurat), levels = c("RPC0","PhRep1","RPC2","RPC3", "NEUR4", "RGC5", 'COR6', "COR7", "NEUR8", "POM9", "RPE10", "AMA11", "RPC12", "POM13", "BV14", "CilEp15", "WBC16", "17", "EOM18", "RBC19", "LSC20", "RGC21", "MEL22", "23", "24", "CilEp25"))
markers.to.plot <- c("Ptf1a","Onecut2","Prdm13","Hbb-y","Hba-a2","Egfl7","Ramp2","Igfbp7","Dct","Mitf","Pmel","Tyr","Pax3","Sox10","Foxd3","Vsx2","Lhx2","Sorbs2","Eno1","Ccnd1","Mgp","Pitx2","Irx1","Klf5","Krt7","Krt5","C1qc","C1qb","Osr1","Klf4","Foxl2","Cdh6","Kera","Matn4","Aqp1","Aldh3a1","Pax7","Myog","Mstn","Ecm1","Rgs5","Tbx18","Mgp","Lhx4","Otx2","Crx","Prdm1","Foxn4","Notch1","Rbpj","Postn","Col3a1","Tnmd","Isl1","Pou4f2","Atoh7")
markers.to.plot <- unique(markers.to.plot) ## needed to fix Error in `levels<-`(`*tmp*`, value = as.character(levels)):factor level [43] is duplicated##
DotPlot(Myrfhumgroup_sub.seurat, features = markers.to.plot, cols = c("blue", "red", "green", "black", "grey", "pink"), dot.scale = 8) +
  RotatedAxis()

Idents(Myrfhumgroup_sub.seurat) <- factor(Idents(Myrfhumgroup_sub.seurat), levels = c("RPC0","PhRep1","RPC2","RPC3", "NEUR4", "RGC5", 'COR6', "COR7", "NEUR8", "POM9", "RPE10", "AMA11", "RPC12", "POM13", "BV14", "CilEp15", "WBC16", "17", "EOM18", "RBC19", "LSC20", "RGC21", "MEL22", "23", "24", "CilEp25"))
markers.to.plot <- c("Ptf1a","Onecut2","Prdm13","Hbb-y","Hba-a2","Egfl7","Ramp2","Igfbp7","Dct","Mitf","Pmel","Tyr","Pax3","Sox10","Foxd3","Vsx2","Lhx2","Sorbs2","Eno1","Ccnd1","Mgp","Pitx2","Irx1","Klf5","Krt7","Krt5","C1qc","C1qb","Osr1","Klf4","Foxl2","Cdh6","Kera","Matn4","Aqp1","Aldh3a1","Pax7","Myog","Mstn","Ecm1","Rgs5","Tbx18","Mgp","Lhx4","Otx2","Crx","Prdm1","Foxn4","Notch1","Rbpj","Postn","Col3a1","Tnmd","Isl1","Pou4f2","Atoh7")
DotPlot(Myrfhumgroup_sub.seurat, features = markers.to.plot, cols = c("blue", "red", "green", "black", "grey", "pink"), dot.scale = 8) +
  RotatedAxis()

##Making dot plot to look at downregulated genes in RPE cells across all clusters
DefaultAssay(Cluster.merged.seurat) <- "RNA"
Idents(Cluster.merged.seurat) <- factor(Idents(Cluster.merged.seurat), levels = c("RPC", "NEUR", "Oligos", "AMA", "RPE", "MEL", "PERI", "CilEp", "WBC", "RBC", "BV", "EOM", "POM", "LSC", "COR", "PhRep", "RGC", "CONJ", "ONH"))
markers.to.plot <- c("Ankfn1",	"Pde1c", "Otx2os1",	"Trpm3",	"Col3a1",	"Serpine3",	"Col2a1",	"Ttn",	"Rora",	"Rbms1",	"Col1a1",	"H19", "Col1a2",	"Sox5",	"Il1rapl1",	"Ccdc141",	"Postn",	"Etl4",	"Dcn",	"Sema3c",	"Zbtb20",	"Cacna1c",	"Kcnq5",	"Mgp",	"Tcf4",	"Auts2",	"Rasgrp3",	"Cdh3",	"Gm45895",	"Pitpnc1",	"Pde4b",	"Nfia",	"Sox9",	"Nfib",	"Col5a2",	"Mir99ahg",	"Eln",	"Ahnak",	"Tshz2",	"Col11a1",	"Ak5",	"Pip5k1b",	"Prdm16",	"Mitf",	"Igf2",	"Gm10076",	"Sntb1",	"Dlgap1", "Lgals1",	"Sox6",	"Gulp1",	"Dlg2",	"Oca2",	"Fndc3b",	"Prkg1",	"Mir100hg",	"Gpc6",	"Fbn2",	"Samd4",	"Slc4a5",	"Neb",	"Airn",	"Maml2",	"Rian",	"Stard13",	"Trabd2b",	"Wfikkn2",	"Rbms3",	"Matn4",	"Filip1l",	"Aff3",	"Lox",	"Creb3l2",	"Kcnq1ot1",	"Bnc2",	"Fbxl7",	"Zfp385b",	"A930028N01Rik",	"Pbx1",	"Uba52",	"Meg3",	"Wwp2",	"Ankrd44",	"Trpm1",	"Dock10",	"Gm48956", "Camk1d",	"Phactr1",	"Ppfibp2",	"Vim",	"Azin2",	"Myh3",	"Myrip",	"Creb5", "Ehbp1",	"Gm26917",	"Pard3b",	"Igf1",	"Capns1",	"Asic2",	"Dock9", "Dlc1",	"Pxdn",	"Peg3",	"Pbx3",	"Fmnl2",	"Magi2",	"Rsph10b",	"Fgfr2",	"Acta1",	"Phldb2",	"Mef2c")
DotPlot(Cluster.merged.seurat, features = markers.to.plot, cols = "RdYlBu", dot.scale = 8) +
  RotatedAxis()

## Making dotplot to look at down DEGs in RPE cells
Idents(Myrfhumgroup_RPE.seurat) <- factor(Idents(Myrfhumgroup_RPE.seurat), levels = c("RPE"))
markers.to.plot <- c("Ankfn1",	"Pde1c", "Otx2os1",	"Trpm3",	"Col3a1",	"Serpine3",	"Col2a1",	"Ttn",	"Rora",	"Rbms1",	"Col1a1",	"H19", "Col1a2",	"Sox5",	"Il1rapl1",	"Ccdc141",	"Postn",	"Etl4",	"Dcn",	"Sema3c",	"Zbtb20",	"Cacna1c",	"Kcnq5",	"Mgp",	"Tcf4",	"Auts2",	"Rasgrp3",	"Cdh3",	"Gm45895",	"Pitpnc1",	"Pde4b",	"Nfia",	"Sox9",	"Nfib",	"Col5a2",	"Mir99ahg",	"Eln",	"Ahnak",	"Tshz2",	"Col11a1",	"Ak5",	"Pip5k1b",	"Prdm16",	"Mitf",	"Igf2",	"Gm10076",	"Sntb1",	"Dlgap1", "Lgals1",	"Sox6",	"Gulp1",	"Dlg2",	"Oca2",	"Fndc3b",	"Prkg1",	"Mir100hg",	"Gpc6",	"Fbn2",	"Samd4",	"Slc4a5",	"Neb",	"Airn",	"Maml2",	"Rian",	"Stard13")
DotPlot(Myrfhumgroup_RPE.seurat, features = markers.to.plot, cols = "RdYlBu", dot.scale = 8, split.by = "Geno") +
  RotatedAxis()

## Making dotplot to look at down DEGs in RPE cells
Idents(Myrfhumgroup_RPE.seurat) <- factor(Idents(Myrfhumgroup_RPE.seurat), levels = c("RPE"))
markers.to.plot <- c("Trabd2b",	"Wfikkn2",	"Rbms3",	"Matn4",	"Filip1l",	"Aff3",	"Lox",	"Creb3l2",	"Kcnq1ot1",	"Bnc2",	"Fbxl7",	"Zfp385b",	"A930028N01Rik",	"Pbx1",	"Uba52",	"Meg3",	"Wwp2",	"Ankrd44",	"Trpm1",	"Dock10",	"Gm48956", "Camk1d",	"Phactr1",	"Ppfibp2",	"Vim",	"Azin2",	"Myh3",	"Myrip",	"Creb5", "Ehbp1",	"Gm26917",	"Pard3b",	"Igf1",	"Capns1",	"Asic2",	"Dock9", "Dlc1",	"Pxdn",	"Peg3",	"Pbx3",	"Fmnl2",	"Magi2",	"Rsph10b",	"Fgfr2",	"Acta1",	"Phldb2",	"Mef2c")
DotPlot(Myrfhumgroup_RPE.seurat, features = markers.to.plot, cols = "RdYlBu", dot.scale = 8, split.by = "Geno") +
  RotatedAxis()


## Validating top downregulated hits in hMYRF RPE cells
DefaultAssay(Myrfhumgroup_RPE.seurat) <- "RNA"
Idents(Myrfhumgroup_RPE.seurat) <- factor(Idents(Myrfhumgroup_RPE.seurat), levels = c("RPE"))
markers.to.plot <- c("Sema3c",	"Dcn",	"Etl4",	"Postn",	"Ccdc141",	"Il1rapl1",	"Sox5",	"Col1a2",	"H19",	"Col1a1",	"Rbms1",	"Rora",	"Ttn",	"Col2a1",	"Serpine3",	"Col3a1",	"Trpm3",	"Otx2os1",	"Pde1c",	"Ankfn1")
DotPlot(Myrfhumgroup_RPE.seurat, features = markers.to.plot, cols = "RdYlBu", dot.scale = 8, split.by = "Geno") + RotatedAxis()

#Making a feature plot to look at PRSS56 in RPC cell cluster
DimPlot(Cluster.merged.seurat, reduction = "umap")
FeaturePlot(Cluster.merged.seurat, features = "Prss56", split.by = "Geno", max.cutoff = 3, cols = c("grey", "red"))
VlnPlot(Cluster.merged.seurat, features = "Prss56", idents = "RPC", split.by = "Geno", pt.size = 0)
VlnPlot(Cluster.merged.seurat, features = "Prss56", idents = "RPC", split.by = "Geno") + ylim(0,4)

#this graph is too busy do RPC cluster only
VlnPlot(Myrfhumgroup_RPC.seurat, features = "Prss56", group.by = "Geno", pt.size = 0.5) + ylim(0,4)
VlnPlot(Myrfhumgroup_RPC.seurat, features = "Prss56", group.by = "Geno")
FeaturePlot(Myrfhumgroup_RPC.seurat, features = "Prss56", split.by = "Geno", max.cutoff = 3, cols = c("grey", "red"))
table(Myrfhumgroup_RPC.seurat@active.ident, Myrfhumgroup_RPC.seurat$Geno)
table(Myrfhumgroup_sub.seurat@active.ident, Myrfhumgroup_sub.seurat$Geno)

#Subset Myrfhumgroup_sub.seurat.Rdata with only hMYRF Genotypes
Myrfhumgroup_hMYRFGenoONLY <- subset(Myrfhumgroup_sub.seurat, subset = Geno == c("Myrf Control", "Myrf Mutant"))
Myrfhumgroup_hMYRF_Control <- subset(Myrfhumgroup_sub.seurat, subset = Geno == ("Myrf Control"))
Myrfhumgroup_hMYRF_Mut <- subset(Myrfhumgroup_sub.seurat, subset = Geno == ("Myrf Mutant"))
plot1 <- DimPlot(Myrfhumgroup_hMYRFGenoONLY, reduction = "umap", label = TRUE) + NoLegend()
plot2 <- DimPlot(Myrfhumgroup_hMYRF_Control, reduction = "umap", label = TRUE) + NoLegend() + ylim(-14,15)
plot3 <- DimPlot(Myrfhumgroup_hMYRF_Mut, reduction = "umap", label = TRUE) + ylim(-14,15)
plot1 + plot2 + plot3
plot4 <- FeaturePlot(Myrfhumgroup_hMYRFGenoONLY, features = "Myrf", split.by = "Geno")
plot5 <- FeaturePlot(Myrfhumgroup_RPE.seurat, features = c("Myrf", "Tmem98", "Ermn"), split.by = "Geno")
plot6 <- FeaturePlot(Myrfhumgroup_hMYRF_Control, features = "Myrf") + ylim(-14,15)
plot7 <- FeaturePlot(Myrfhumgroup_hMYRF_Mut, features = "Myrf") + ylim(-14,15)
VlnPlot(Myrfhumgroup_RPE.seurat, features = c("Myrf", "Tmem98", "Ermn"), group.by = "Geno", cols = c("Grey", "Magenta")) + ylim(0,2.5)
VlnPlot(Myrfhumgroup_RPE.seurat, features = "Ermn", split.by = "Geno") + ylim(0,2.5)
#Saving hMYRF Geno subset
save(Myrfhumgroup_hMYRFGenoONLY, file = "Myrfhumgroup_hMYRFGenoONLY.seurat.Rdata")

#Subset Myrfhumgroup with humMUT and P0 KO timepoint
Myrfhumgroup_hMYRF_P0_Geno <- subset(Myrfhumgroup_sub.seurat, subset = Geno == c("Myrf Control", "Myrf Mutant", "P0_WT", "P0-MUT"))
save(Myrfhumgroup_hMYRF_P0_Geno, file = "Myrfhumgroup_hMYRF_P0_Geno.seurat.Rdata")